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Regularized Sampling of Multiband Signals 

J. Selva 



Abstract — This paper presents a regularized sampling method 
for multiband signals, that makes it possible to approach the 
Landau limit, while keeping the sensitivity to noise at a low 
level. The method is based on band-limited windowing, followed 
by trigonometric approximation in consecutive time intervals. 
The key point is that the trigonometric approximation "inherits" 
the multiband property, that is, its coefficients are formed by 
bursts of non-zero elements corresponding to the multiband 
components. It is shown that this method can be well combined 
with the recently proposed synchronous multi-rate sampling 
(SMRS) scheme, given that the resulting linear system is sparse 
and formed by ones and zeroes. The proposed method allows one 
to trade sampling efficiency for noise sensitivity, and is specially 
well suited for bounded signals with unbounded energy like those 
in communications, navigation, audio systems, etc. Besides, it 
is also applicable to finite energy signals and periodic band- 
limited signals (trigonometric polynomials). The paper includes 
a subspace method for blindly estimating the support of the 

\ multiband signal as well as its components, and the results are 

, validated through several numerical examples. 

i I. Introduction 

Sampling is the operation that makes the discrete processing 
of continuous signals possible. The basic tool for this operation 
is Shannon's Sampling Theorem, which states that a signal can 

■ be reconstructed without error from its samples taken at a rate 
' equal to twice its maximum frequency (Nyquist rate), [1]. In 

some situations, however, the signal is multiband, in the sense 
I that its spectral support is composed of a finite set of disjoint 

■ intervals. In this case, sampling at the Nyquist rate can be 
very inefficient. The design of alternative sampling schemes 

■ for this kind of signals has been an active research topic for 
I decades. It was first demonstrated by H. J. Landau in [2] that if 

■ the signal is sampled nonuniformly, then the average sampling 
rate is lower bounded by the measure of its spectral support, 

1 (Landau lower bound). Later, it was shown in [3] that there 
exist nonuniform sampling patterns that approach this bound 
as much as desired. The sampling scheme in this reference was 
the so-called "multi-coset sampling", which consists in select- 
ing a sampling rate above the Nyquist rate, and then choosing 
a periodic nonuniform subsequence of it. The subsequent 
investigations on this topic have mainly focused on devising 
low complexity and stable (low sensitivity) implementations of 
multi-coset sampling, [4]-[7]. Recently, compressed sensing 
techniques have been applied to multi-coset sampling in [8], 
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[9], so as to detect the band positions, as well as to reconstruct 
each of the components. An alternative method to discretize 
multiband signals is the random demodulator in [10], [11], in 
which the complexity of the analog processing is increased, 
with the aim of employing only analog-to-digital converters 
of moderate input bandwidth. Finally, it is worth mentioning 
the synchronous multi-rate sampling (SMRS) scheme in [12]. 
Here, the multiband signal is sampled in several uniform grids 
with different rates, in order to produce a bunch of lowpass 
signals, in which the spectrum of the multiband signal appears 
aliased. The key idea is that this aliasing may not affect all 
frequencies of the original spectral support for all sampling 
rates. This fact is exploited in [12] by properly selecting 
the number of sampling frequencies, so as to achieve the 
reconstruction of the multiband signal. This technique has been 
further developed in a second paper, which also includes the 
application of compressed sensing [13]. 

An assumption that is common to these approaches is that 
of finite energy: the multiband components are assumed to 
be in the Lebesgue space L^, and the machinery of Fourier 
analysis for this space is applied extensively. This way to 
proceed simplifies the analysis, because signals in have 
a proper spectrum function, and the well-known properties of 
the Fourier transform apply. Nevertheless, there is a mismatch 
between this model and the situation encountered in multitude 
applications, (communications, audio, navigation, etc). Quite 
often, it is not possible to sample a time interval that contains 
most of the signal's energy due to its long duration. This 
limitation has negative effects, that can be readily noticed 
in the single-band case by analyzing the Sampling Theorem. 
Given a finite energy signal s(t) whose spectrum lies in 
[— 1/(2t), 1/(2t)], the Sampling Theorem provides the rep- 
resentation 

oo 

s(nr + u) = s((n — p)r)sinc(p + u/t), (1) 

P— — 00 

where n is an integer and u is any time shift following — t/2 < 
u < t/2. Yet in practice, the infinite set of samples s((n— p)r) 
is not available, and then ([T]| must be truncated at some index 

P 

p 

s(nr + u) w s((n — p)T)sinc(p + w/r). (2) 

p=-p 

This last formula is well known for its poor accuracy, [14]. 
The problem lies in that ([T]i converge only because the samples 
s((n — p)r) converge to zero as \p\ oo, due to the fact that 
s{t) has finite energy. So, if the energy of s{t) is mostly outside 
the sampling interval [—PT, PT], the accuracy is poor due to 
the sine tails corresponding to the neglected samples. 

This problem is well known in the signal processing field, 
and it can be solved by regularizing the interpolation formula 
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in (|2]i, using various filter design techniques, [15]. A typi- 
cal method consists in assuming that there is some excess 
sampling bandwidth, i.e, that the spectrum of s{t) lies in 
[-B/2,B/2] with Bt < 1 (and not Bt < 1), and then 
multiplying the sine samples in (|2]) by a set of weights Vp(u), 

p 

s(nT + u) w s{{n ~ p)T)vp{u)smc{p + u/t). (3) 

p=-p 

This regularization converts the ill-behaved interpolator in (|2]i 
into a well-behaved one, and this can be readily seen in the 
error trends: while the interpolation error of ^ decreases only 
as 1/P, [14], there are sets of coefficients Vp{u) for which the 
error of ^ decreases exponentially as e"'^'^^^'^^^, [20]. The 
price paid for the regularization is that it is not possible to 
represent signals using the full sampling bandwidth anymore, 
and it is also not possible to implement filters with sharp 
transitions, since this would produce convergence problems 
similar to those of the sine series. 

By analogy with the single band case, the multiband prob- 
lem is currently in the stage equivalent to that of the sine 
series in ([T]i and its truncated version in (|2|i; i.e, some sam- 
pling schemes have been identified that approach the Landau 
rate, and allow the recovery of the multiband signal and its 
components, but there is no regularization procedure available. 
This lack of regularization can be noticed in analytical devices, 
like the slicing of the spectrum using sharp transitions, and in 
the use of discrete sequences with full spectrum, like coset 
sequences. 

The purpose of this paper is to present a regularized sam- 
pling scheme for multiband signals. Relative to the existing 
methods in the literature, the one in this paper provides a 
flexible interface between the multiband signal and its discrete 
representation, that allows one to trade sampling efficiency for 
noise sensitivity. In the scheme, the reconstruction procedure is 
specified in the form of weighted trigonometric polynomials, 
which are able to interpolate the multiband signal and its 
components in consecutive time intervals. 

The starting point in the derivation of the scheme is the 
boundness assumption which substitutes the usual finite energy 
assumption, i.e, the multiband components are assumed to 
have bounded time amplitude, but their energy can be unlim- 
ited. To make no assumption about the energy is meaningful 
since, for long signals, it usually lies outside the sampling 
interval and can be arbitrarily large. And to work assuming 
bounded components reflects clear physical and engineering 
constraints: signals have bounded peak power due to restric- 
tions like the maximum transmitted power, or due to systems 
like automatic gain controls. With the new assumption, the 
multiband components have a lowpass equivalent whose real 
and imaginary parts belong to the Bernstein space B^g, [16, 
chapter 6], i. e, they can be described as entire functions of 
exponential type that are real and bounded on the real axis. 
This new description has the drawback that the basic tool 
in spectral analysis, namely the Fourier transform (Plancherel 
theorem), is not directly usable. So, basic analytic procedures 
in the literature on multiband sampling [3], [10]-[13], [17], 
like the division of the spectrum in disjoint intervals, the use of 



discrete sequences with full spectrum (like coset sequences), 
and the use of lowpass filters with sharp transitions are not 
valid on bounded band-limited signals. However, these tools 
can be substituted by a simple yet powerful analytical device, 
as shown in this paper. If z{t) is the multiband signal, the 
device consists in multiplying z{t) by a bandlimited window 
w(t), that approximately "selects" a given finite interval, in 
the sense that its time-domain content is concentrated 
in it. The key result is that there exist windows w{t), for 
which the product z{t)w{t) is a trigonometric polynomial with 
negligible error in the interval selected by w{t). Besides, the 
trigonometric polynomial "inherits" the multiband structure 
of z{t), i.e, its coefficients are formed by bursts of values 
corresponding to the bands of z{t). Sec. |lll] is dedicated to 
this analytical device. 

After this, the design of a sampling scheme is addressed 
in Sec. IIVI where it is designed for sparse trigonometric 
polynomials, since the product z(<)w(t) can be viewed as a 
polynomial of this type with negligible error. In this section, 
it is shown that a finite version of the SMRS scheme produces 
a sparse linear system, in which the sensitivity to noise can 
be reduced by increasing the number of samples. Afterward, 
Sec. |V] shows how the finite sampling grid for z{t)w{t) can 
exactly match an infinite SMRS scheme for z(<). This makes 
it possible to interpolate z{t) and its components at any t. 
The interpolation error is then analyzed in Sec. |VI] Finally, 
the blind sampling problem is studied in Sec. IVIII where 
the MUSIC algorithm from direction of arrival estimation is 
adopted. 

Since the notation in the paper is extensive, the basic 
symbols and operators are described in the next sub-section, 
and there is a list of symbols at the end of the paper in Table 
H] The next section sketches the sampling method, in order to 
give a broad view of the concepts involved. The novel aspects 
of the paper are, fundamentally, the contents of Sees. Hill to [Vl 
and Sec. IVlIl 

A. Notation 

Several conventions have been adopted in the paper to 
simplify the notation: 

• Definitions are performed using the operator "=". 

• R{t, T) denotes the interval 

R[t,T) = [t-T/2,t + T/2[ (4) 

for arbitrary r and T > 0. 

• The symbols J, J', Izw and Isw.m denote sets of distinct 
indices (integers). In I^w and Isw.m the subscripts 'zw' 
and 'sw' remind one of the meaning of these sets: thus, 
Izw and Isw,m are the index sets of z{t) and Sm{t), 
respectively, after windowing using w{t). 

« I J| denotes the number of elements of a set J. 
> The operator " " " denotes periodization with period T, 
that is, for a signal s(t), it is 

oo 

s(t)= J2 <t+pT)- (5) 

p— — OO 
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• Vectors and matrices are written in bold font and in lower 
and upper case, respectively, (m, M). 

• I denotes the identity matrix. 

• [A] /i n denotes the h, n element of matrix A; [p]r denotes 
the r-th element of vector p; and [A]. „ the n-th column 
of a matrix A. 

II. Sketch of the proposed method 

Consider a multiband signal z(t) formed by M bounded 
band-limited components Sm(i), 

M 

z{t) = «™(*)' 

m— 1 

where the band of Sm{t) is [am,bm], and these bands are 
known and appear in increasing order, bm < a„i+i, m = 
1, . . . , Af — 1. The m-th component can be viewed as a signal 
whose baseband version 

s,n{t)e~^^^-^+'>'-^' (7) 

has real and imaginary parts in the space _^ y The 

spectral support of z{t) is the set 

M 

Sz= [J [am.bm]- (8) 

m— 1 

The multiband sampling problem can be posed in a finite time 
interval, centered at an arbitrary instant r, R{t,T), T > 0, in 
which a set of samples is taken at distinct instants r + ii, r + 
t2, ■ ■ ■ ,T + tpf, but assuming that the components Sm{t) must 
be interpolated only inside a shorter interval R{t, Ti), Ti < T. 
The objective is then to find a set of functions h„j.„(t), such 
that the components Sm{t) and z{t) can be interpolated using 

N 

Sm(T + t) w ^z(r + t„)h™,„(t) (9) 

n=l 

and 

M N 

z{T + t)^J2J2 <^ + (10) 

m— 1 n— 1 

for t in i?(0,Ti). If there is a method to solve this last 
problem, then it can be applied repeatedly in overlapping 
intervals R{t + kTi,T) with integer fc, so as to interpolate 
z{t) and the components Sm(t) at any t. So, for example, the 
value of any of the components Sm{t) in a regular grid of 
instants nAt, with integer n and arbitrary At > 0, could be 
obtained by evaluating (|9|l in the intersection of this grid with 
R{t + kTi,T), and repeating the procedure for consecutive 
values of k. 

In this setting, consider a band-limited window w{t) with 
spectrum lying in [—Bw/2, Bw/2], Bw > 0, that is approxi- 
mately time limited to the interval R{0,T). Besides, assume 
that it is bounded between one and a minimum value greater 
than zero in i?(0,ri), that is, 5^ < w(t) < 1 in i?(0,ri), 
5^ > 0. Note that w{t) is not the typical window in filter 
design, that sharply selects a given time interval and nulls the 
rest of the time axis. Actually, it is the dual case: the spectrum 
of w{t) sharply selects the band [—Bw/2, Bw/2], and w{t) 



is very small but not strictly zero outside R{0,T). If B^ is 
smaller than the minimum separation among bands of z{t). 

By, < mmarn+i - bm, m = 1, . . . , Af - 1, (11) 

then the product z(t + t)w{t) is also a multiband signal with 
components Sm(r + t)w{t), 

M 

z(r + i)w(t) = ^ s„(T + i)w(i), (12) 

m—l 

and spectral support 

M 

5,^ = IJ [a„-B„/2,6„ + B„/2]. (13) 

m— 1 

Thus, z(t + t)w{t) is another multiband signal in which 
the bands of z{t) have been expanded, but are still disjoint 
due to (fTTT i. Now, the multiband sampling problem can be 
posed in (fT2l l with sampling instants i„ and sample values 
z(t + tn)w{tn)- If there is a satisfactory solution for this 
last problem, then the original multiband signal z{t) and its 
components s„i{t) can be readily interpolated in R{t,Ti), 
simply by dividing by w{t). 

In the windowed model in (fT2] i. the signal z(t + t)w{t) 
is approximately time-limited. This allows one to apply the 
Sampling Theorem, but with the time and frequency domains 
switched, i.e, z{t + t)w{t) is repeated periodically with period 
T in the time domain, and this operation samples the fre- 
quency domain. Here, the window w{t) acts as a time-domain 
"lowpass filter" with nonuniform response. The result of this 
operation is a trigonometric polynomial, with coefficients 
equal to samples of the spectrum of z(t + t)w{t), taken at 
the frequencies p/T with integer p, that lie in (fTSl l. If Isw,m 
denotes the integer indices p of the frequencies p/T inside the 
band of Sm(T + t)w{t), 

h^,m = {p : am-B^/2 < p/T < b^+B^/2,p e Z}, (14) 

and Izw denotes the union of these sets, 

AI 
m—l 

then this approximate sampling theorem says that z{T + t)w{t) 
is the polynomial 

z(r + i)w(i)« Cp(r)e-'"2-P*/^, (16) 

for t in i?(0, T), where the Cp(r) are unknown samples of the 
spectrum of z{t + t)w{t). Besides, the components Sm{t) can 
be interpolated using 

pet:,™ 

for t in i?(0,Ti). 

As will be shown in this paper, the accuracy of this approx- 
imation depends on the window w(t). Actually, it is shown 
in the next section that the interpolation error is bounded by 
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Fig. 1. BPSK signal (bandwidth 136/T, roll off 0.8, chip rate 0.01324T), 
window w(t) (B^T = 13.6), and windowed signal (bold curve) in Sec. 
IVIII-AI 



polynomial. This implies that the windowed problem in Eq. 
( fT2b can be cast in a generic setting, as that of computing the 
coefficients /3p of a sparse trigonometric polynomial from its 
samples at the abscissas ti,t2, ■ ■ ■ ,tN. If J is an index set, 
then this kind of polynomial can be expressed as 

The polynomial on the right in ( fTSI l has this form with J — I^w 
and f3p = Cp(r). In this generic setting, the sampling problem 
reduces to solving the linear system 

a(i„; J) = ^/?pe^2.pt„/r^ n=l,...,7V. (20) 

peJ 

If this system has full column rank, then the coefficients /3p 
can be obtained using 



N 



(21) 
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Fig. 2. Interpolation en'or in validation of the interpolation method in Sec. 

ivni-Al 



where 77^ „ denotes the pseudo-inverse of ( |20] l. Then, the 
solution of the sampling problem for a{t: J) is obtained by 
substituting this last equation into ( fT9l ). This yields 



N 



ait-, J) ^J2a{t^;J)e^{t;J) 



n=l 



where 



(22) 



(23) 



Besides, since any "component" of a{t; J) is specified by a 
subset J', J' C J, its associated sampling formula is obtained 
by replacing 6n{t; J) with 6'„(i; J') in ( l22b . 



N 



a{t;J') = ^a(i„; J)0„(t; J')- 



(24) 



yl^e/iJ^u,, where is a bound on the amplitude of z{t), and 
e a bound on the summation 



J2 |w(t + pT)| <e, Uni?(0,T). 

|p|#o 



(18) 



The fact is that there exist windows, as the one presented 
in this paper, for which S^j grows with T, and e decreases 
exponentially as T — > 00 with trend roughly equal to c^^^™^. 
Thus, there is a moderate value of T for fixed such that 
([TTT i can be regarded as an equality, for any fixed numerical 
precision. There are several methods to construct this kind 
of window, [18]-[20]. The one proposed in this paper in Sec. 
IIII-Al is based on the properties of the Fourier transform of the 
Kaiser-Bessel function. Fig. [T] shows this kind of interpolation 
for a baseband BPSK signal s(t) using the window in this 
paper. The bold curve is the windowed signal, that is selected 
for e — 10^* in ( fTsl l. Fig. |2] shows the interpolation error, 
which is below -200 dB roughly in i?(0,r/2). (For details, 
see Sec. IVIII-AI ) 

Coming back to the formula in ( fTSI l, if e is sufficiently 
small, then z(r + t)w(t) can be regarded as a trigonometric 



Finally, the solution for a sparse polynomial in jTzi can be 
applied to the windowed model in ( fT2] l. if J' is identified 
with each of the sets Isw.m, and J with the full set of 
spectral samples I^w So, if (l22b is used with sample values 
z{t + tn)w{tn) and index set Isw.m, and then the effect of the 
window in (fT2l l is removed by dividing by w{t), the result is 
an interpolation formula for Sm (t) 

1 ^ 

(r + i) « — — > z(t + tn)w{t 

n)^ nit ^ I sw .m) : (25) 

for t in -R(0, T\). And the same can be done for interpolating 
z(i), but with index set l^w, 

1 ^ 

z(t + <)«——> z(r + t„)w(i it-hw). (26) 
w(<) ^ 

^ ^ n—l 

The last two formulas solve the initial problem in (|9]l, with 
functions „(t) given by 



(t) =witn)0nit;ls )/w(0. 



(27) 



This is in broad terms the sampling method proposed in this 
paper for multiband signals, and the next three sections are 
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dedicated to clarify its fundamental aspects. The windowing 
and spectral sampling are explained in the next section. The 
key points in it are the approximate truncation of the multiband 
signal using w{t), and the dual sampling theorem, which is 
demonstrated by means of the Poisson's summation formula. 
Also, a specific window w{t) is selected in Sub-section IIII- Al 
Then, Sec. |IV] studies the selection of the instants t„ so that 
the linear system in ( |20l l has low sensitivity to perturbations. 
It turns out that a finite SMRS scheme produces a sparse 
linear system of ones and zeroes, in which it is possible to 
reduce the sensitivity by slightly over-sampling. Next, Sec. 
IVIII shows that the finite SMRS scheme in Sec. |IV] can be 
perfectly integrated into an infinite SMRS scheme for the 
initial multiband signal. This means that the finite scheme in 
an interval R{t, T) with arbitrary r employs samples lying in 
the infinite scheme exclusively. 

III. An approximate dual sampling theorem: 

APPROXIMATION BY MEANS OF TRIGONOMETRIC 
POLYNOMIALS 

Recall the multiband signal z(t) in ^ with components 
Sm(i), and let Ag^m denote specific bounds for them. 



and 



(28) 



Also assume that it is necessary to approximate the value of 
Sm{t) for some or all m, 1 < m < M, in an interval R{t, Ti) 
for a specific r. To perform this approximation, take a band- 
limited window w(i) in with spectral support contained 
in [-B„/2, i?^/2], that is bracketed in i?(0,ri) as 5^ < 
w{t) < 1 for a specific 6^ > 0, and that additionally fulfills 
the concentration property in Eq. (fTsT l with Ti < T. Next, 
form the new signal 



,{t;T) = z{T + t)w{t). 



(29) 



This product expands each of the bands [am, bm] of z{t + t), 
so that the spectral support of Zw{t; t) is contained in the set 
( fTsT l. Besides, this windowing affects all components s,„(t) 
of z(t) equally, i.e, z^(i;T) also consists of M components 

M 

Zw{t;T) = ^ s^,m(i;r), (30) 



where 



Stu,m(i; t) = Sm{t + r)w(t). 



(31) 



If is smaller than the minimum separation among bands 
[Eq. (fTTl il. it turns out that z^{t; r) is also a multiband signal 
(its bands are disjoint). The inequalities in (fTsl l and ( |28] | imply 
that Sui,m{t; t) and z^{t; r) can be well approximated by their 
respective periodic versions, Zu,(i;r) and Su),m(i;T), defined 
using (|5]). It is 



= I 51 + * + pr)w(t +pT)\< As^^e 

p#0 



(32) 



M 



\z.u,{t;T)~z{t)\ < ^ |s.u,,„i(<;r)-Su,,m(i;T)| < e ^ ^s,m- 

m— 1 m— 1 

(33) 

Now, the signals Su,.m(i;r) and Zu,(i;T) are band-limited 
and periodic. So, by the Poisson's summation formula [16, 
Sec. 2.3], they are sparse trigonometric polynomials, whose 
coefficients correspond to nonzero samples of the spectra 
of Su),m(i; t) and z,„(t;T) at frequencies of the form p/T, 
respectively. The sets of indices p of these frequencies, denoted 
Isw,m and Izw, have akeady been defined in (fT4l) and (fTSl) . 
Thus, Sw,m{t\T) and z^(i;T) can be expressed as 



and 



(34) 



(35) 



where Cp(r) is the (unknown) value of the spectrum of 
z^(t;r) at frequency p/T. So, the disjoint spectral intervals 
[am,bm] of z{t) translate into disjoint bursts of coefficients 
in Zuj(t;T). Note that Su,,m(t;T) and z^(t;r) have proper 
spectrum functions, because they are band-limited signals that 
lie in L^. This is not true for either Sm{t) or z(t). 

Finally, the desired formulas are obtained from ( |34] | and 
For t in i?(0,ri), it is 



j2Tvpt/T 



with error bounded by Ag^me/Sw, and 

+ « — = — ^ y] Cp(r)i 



(36) 



with error bounded by {J2,n ^s,m)^/Sw 

A specific window w{t) is presented in the next sub-section. 

A. Selection of a window function 

As shown in several apphcations [20]-[23], a band-limited 
window with excellent concentration properties can be con- 
structed from the Fourier transform of the Kaiser-Bessel func- 
tion. For the interpolation problem in this paper, the proposed 
window is 



^^^^ ^ sinc(J5^0sinc((l - S)B^^t^ - p^T^/A) ^^^^ 



sinc(j(l ~ 5)pB^T/2) 



where 



p = v/1 - i/(Bu.r)2 



(39) 



In ( 138b . the second sine function in the numerator is the Fourier 
transform of the Kaiser-Bessel function, and the concen- 
tration of w{t) is roughly proportional to the denominator in 
(|38] l, which increases exponentially as e'^'^^*'''^'"-^/^/r when 
T — oo. The first sine function in ( [38] l is necessary so as 
to damp the tails of w(i), given that it would otherwise not 
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Fig. 3. Error bound t/Sw using the values for 5 and e in )4Qt , for Ti 

3T/4,T/2,T/4. 



belong to L^. A bound e for the inequality in (fTST l is derived 
in Ap. H] for this window, which is well approximated by 



^opt ' 



! 0.03326 - 0.002084B^T + 0.3737 • 10"'*(B,„r)^ 



il. 086-0. 6676B„T 



(40) 

Here, 5opt is the value of 6 in ( l38T l that achieves the cor- 
responding e in ( |40] |. Note the exponential rate in this last 
expression for e. As a consequence, the model accuracy can 
be increased arbitrarily by increasing B^T. For B^T = 13.61 
it is e = 10"** and for B^T = 25.59 it is e = lO'^^. The 
actual bound on the interpolation error is given by e/S^, where 
(5tu depends on Ti . A value for 5^ can be obtained numerically 
using the definition 



inf w(i). 

|t|<Ti/2 



(41) 



The error bound e/S^ versus B^T is plotted in Fig. |3] for 
Ti = 3r/4,T/2,r/4. In this figure 5 and e were selected 
using (l40l i. Note that any practical accuracy can be achieved 
by increasing the product B^T. 

IV. Low SENSITIVITY SAMPLING OF A SPARSE 
TRIGONOMETRIC POLYNOMIAL 

The generic sampling problem for sparse trigonometric 
polynomials, as posed in Sec. Eqs. (fT9] l and ( |20] |. does 
not impose any constraint on the distribution of the instants 
i„. However, the condition number of ( |20l i can vary wildly 
with it, as can be easily checked numerically. The distribution 
in which the sampling points are equally spaced and cover 
i?(0, T) is interesting, because it exploits the sparseness of 
the polynomial, i.e, it produces equations in which most 
coefficients are zero. 

To see this point, assume a{t\ J) is uniformly sampled at 



tq+l — to 



It 



(42) 



for q = 0, 1, . . . , Q — 1, Q > 0, where to is an arbitrary instant, 
and define the coefficients 



Next, compute the DFT of the sequence a{tq^i;J) (divided 
by Q), defining a corresponding set of coefficients A^: 

-, Q-i 
, Q-1 



9=0 pGJ 
1 Q-1 

^ 9=0 pSJ 



p{to/T+q/Q)^-j2^rq/Q ^^^^ 



Q-1 



^ 9=0 



{p~r)q/Q 



peJ 

p=r (mod Q) 



peJ 

p=r (mod Q) 



This derivation shows that A^ can be computed either from 
the samples a{tq+i; J) using its definition. 



1 ^ 



-j2TTrq/Q 



or from the coefficients 6,p using 



Ar 



peJ 

=r (mod Q) 



(45) 



(46) 



So, if the samples Q!(iq_|_i; J) are known, Eq. ( |45] | allows 
one to compute the A^ with the best possible conditioning, 
given that DFT matrices have condition number equal to one. 
And Eq. (|46] | shows that A^ is just the sum of some of the 
coefficients 5p (those with index p congruent with r modulo 
Q). In some cases, the sum in Eq. ( l46b may contain a single 
element, and then one of the coefficients Sr would be revealed 
by Ar, i.e, 6p = A,, with p = r(mod Q) and < r < Q. 
Usually, however there are several summands, but if the A^ 
are known, then (l46T l is in any case a sparse linear system of 
ones and zeroes. 

This argument suggests that a regular grid like ( |42] | is a 
good choice, in order to obtain a well conditioned linear 
system. But the selection of Q is problematic. On the one 
hand, if Q is larger than or equal to the difference between 
the largest and the smallest element in J, then the sum in ( l46t 
has a single summand, and the A^ reveal the coefficients Sp 
completely. But the over-sampling would be too large in this 
case. On the other hand, if Q is smaller than that difference, 
the system in ( l46l l could be under-determined. The solution 
proposed in the sequel consists in adopting a finite version of 
the SMRS scheme in [13], in which several integer moduli 
Qi,Q2, ■ ■ ■ ,Qk are used instead of a single one. 

Let us state first the solution of the interpolation problem 
for this scheme, and then study its sensitivity. For a fixed set 
of moduli Qk, let t^^q+i denote a finite SMRS scheme with 
initial instant to, 



(43) 



tk,q+i =to + Tq/Qk, g = 0, 1, . . . ,Qfe - 1, 



(47) 
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and let denote the scaled DFT for modulo Qk as in ( l44l i. 

Then, Eq. (|46] | for all moduli Qfe can be written as 

Afc,. - ^ (49) 

p=r (mod Qfc) 

fc = \,2, . . . ,K, q = 0, 1, ... , Qfc. This is a linear system 
with unknowns Sp. If it has full column rank, then its solution 
can be written is terms of its pseudo-inverse, denoted \^k,r, 
i.e, 

K Qfc-l 

5p — ^ ^ \,k,r^k,r (50) 
k=l r=0 

Now, Eqs. (l43T l and (fT9] l allow one to write a(i; J) in terms 
of the coefficients 6p, 

a(t;J) = ^V'2"^'*"*°^^^- (51) 

Finally, if Eqs. (l48T l and dSOl l are substituted into this equation, 
and the summations are then reordered, the result is 

K Qk-1 

ait; J) = ^ ^ a{tk,g+i;J)9k,g{t;J,to), (52) 

k=l q=0 

where 

(^k,q{t; J, to) 

. Qk-1 

pGJ r— 

(53) 

The last two equations specify the solution to the interpolation 
problem using the finite SMRS scheme. Note that ( l52l l can 
be used to reconstruct any component of a{t; J), simply by 
restricting the index set in 9k.q{t\ J, to), i.e, if J' C J, then 

a{t;J')^J2 J2 a{tk,q+i;J)ek,g{t;J',to). (54) 

fc=l q=0 

Also, note that here 9k.q{t; J,to) is the specific reconstruction 
function in the SMRS scheme with initial instant to, while 
On{t; J) in Sec. |II] Eq. (l22T l. is a reconstruction function for a 
generic sampling set ti,t2, ■ ■ ■ ,tN and fixed r. 

The selection of the moduli Qk can be done now so as to 
make this method viable, i.e, so that the linear system in (|49] l 
has full column rank, and the sensitivity to perturbations of 
(I52I 1 is low. Let us define first the sensitivity measure 

K Qk-1 

7(i; J', to) = (E E to)\y^'- (55) 

k=l q=a 

and its supremum in R{0, T), 

7(J')= sup 7(<;^',0). (56) 

-T/2<t<T/2 

7(i; J', to) is the standard deviation factor due to any pertur- 
bations in the samples a{tk^q\J')- So, if these samples are 



contaminated by perturbations of zero mean and variance a^, 
then the value of a{t\ J') computed using (|54| | has variance 

aV(^;^',^o). 

The following method gives a suitable collection of moduli 

Qk- 

1) Select positive integers Qi and K following QiK \ J\ 
and then sti Qk — Qi + k — 1, k — 2, . . . . 

2) Check whether the linear system 

p=r (mod Qk) 

is over-determined (or determined). If not increase either 
K or Qi. 

3) Finally, the sensitivity can be reduced by adding other 
Q which are selected so as to minimize 7(J), or 7( J') 
if the performance for a specific components J' C J 
must be improved. 

V. Interpolation of a bounded multiband signal 

USING THE SMRS SCHEME 

As shown in Sec. HUl Eq. (l37l i. the sampling problem for a 
bounded multiband signal z(t) can be reduced to a problem 
of the same type, but involving a sparse trigonometric poly- 
nomial. Specifically, it was shown that there is a polynomial 
z(t; r) such that 

+ (58) 

where t lies in R{0,T). This formula specifies how the 
multiband signal z(<) can be interpolated if the polynomial 
z{t:T) is known, and how the polynomial z{t;T) can be 
obtained from samples of z{t). 

The objective of this section is to show that the finite SMRS 
scheme in the previous section of the form 

to+Tq/Qk, g-0,l,...,Qfc-l, k = l,2,...,K, (59) 

can be integrated into the infinite SMRS scheme 

nT + Tq/Qk, n in Z, g = 0, 1, . . . , - 1. (60) 

More precisely, there is a to such that there is a one-to-one 
correspondence between the samples of z{t) in the scheme 
(I6OI 1 that lie in R{t,T), and those of z(t;r) in the scheme 
( |59] l. To demonstrate this, assume that it is necessary to sample 
z{t;T) at one of the instants in ( |59] l. Let nk,q{to) denote 
the unique integer such that to + TqjQk + ^k q{to)T lies in 
i?(0,T), 

nk,qito) = \iT ~ to)/T - 1/2 - q/Qk^ . (61) 
From ( ISST l and recalling that z{t; r) is T-periodic, it is 

z{to + Tq/Qk-.r) = z{to + Tq/Qk + nk,q{to)T;T) 
^z{T + to+Tq/Qk + nk,q{to)T) (62) 
• w{to + Tq/Qk + nk^q{to)T) 



g 



The argument of z in this formula belongs to the scheme in 
(|60] l for any pair of indices (fc, q) if — —t. So, for this 
it is 



z(-T + Tq/Qk\ t) « 7.{Tq/Qk + nk,q{-T)T) 

■ w{~T + Tq/Qk + nk,g{-T)T). 



(63) 



This is the one-to-one correspondence between the samples 
of z{t; t) in ( |59] l and those of z(t) in the intersection of ( |60] l 
with R{t,T). 

Now, it is possible to write an interpolation formula for 
z{t). For this it is enough to "sample" z{t) using (l63T l. then 
obtain the polynomial z{t) approximately using the formula in 
(|54| |. and finally undoing the windowing [dividing by w{t)]. 
To simplify the notation, define first the instants 

tlgir) = -T + Tq/Qk + nk,g{-T)T. (64) 

From (|63]) . the samples of z{t\ r) are 



~z{-T + Tg/g^; r) « z(t + i'fe,,(T))w(i'fe, Jr)). 



(65) 



If the sampling formula in ( l54l i is applied to z(t; r) with these 
sample values, it is 

K Qfc-l 

XI z(T + 4_,(r))w(4_^(r))6'fc,9(t;/^™,-r) 

(66) 

And, finally, ( |58] ) gives the desired formula for z(f), 

K Qfc-1 



^ ^ fc=l g=0 



(67) 



The formula for s„i{t) is the same but with replaced with 



VI. Error analysis 



In the error analysis in the sequel, the sampling scheme is 
denoted using a single index as ti,t2, ■ ■ ■ ,tN instead of two, 
so as to simplify the notation. This also affects the function 
9k,q{t, J, to) in dSSl l, that will be denoted 6'„(t, J, to). Also, the 
error terms will be written inside square braces for readability. 

Let us recall the sampling and interpolation process on z{t) 
in the last three sections, but keeping track of the different 
errors, so as to produce an error bound. First, the multiband 
signal z{t) can be approximated using the polynomial z^(t; r), 
[Eq. mi 



-t) 



1 



w(t) 



iw{t; t) + 



1 



w(t) 



{zyj{t;T)~Zw{t;T)) 



(68) 



Next, Zyj {t; r) can be reconstructed using the sampling formula 
in (|52] | with to = — t, J ~ Iziv, and samples Zw{tn] t), 



N 



(69) 



The samples Zu,(i„;T) are approximately equal to z(t 



Zw{tn; t) = z(t + t„)w(t„) + z^(t„; t) - Zw{tn; t) 



(70) 



But z(T + t„) may be contaminated by a perturbation, denoted 



z{t + tn) = z{t + tn) + ri{tn) + - r]{t„) 



(71) 



Next, replace in dTOl l the term z(t + 1„) with the right side of 
(ItTT i. so as to obtain 



Zu,(t„; r) = (z(t + tn) + ri(tn))w(tn) 



- 1l{in)ff{tn) + iw{tn] t) - Z^(t„; t) 



(72) 



In turn, substitute this equation into (I69l l. and finally use the 
result to replace the term z^(i;r) in (|68] |. The outcome of 
these replacements is a formula with error term for z(t + i): 



1 ^ 

!(r + <) = — r + i„) + r7(t„))w(i„)6'„(t; hn,.~T) 

' n=l 



+ 



N 



-* n=l 

H — 7;T(z«,(i;r) - z^(t;T)) 
w(t) 

(73) 

Now, since w(t) < 1 in i?(0, Ti), the three error terms can be 
bounded in the following way. For the first term in ( |73]) : 



N 



^ ' n—1 ^ 

< — 7777(^;4,™,-t), 



(74) 



where 



and 



n=l 



1/2 



1/2 



(75) 



(76) 



For the second. 



n=l 



(77) 



w(t) 



m—l 



And for the third, 
1 



M 



w{t) 



(z^(t;r) - z^,(t;r)) 



< 



^ ' m—l 



E^^ 



(78) 
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So, the final bound is 

z{T + t) 



N 



< 



w{t) ^ 

Ar,^(t\Iy_^y,,-T) 



w{t) 



M 



(79) 

The first term is due to the sensitivity of the linear system. 
The second is due to the mismatch and perturbations, and it 
can be arbitrarily small, since e tends to zero exponentially 
with T, [Eq. 

The error analysis can be repeated for the components 
Sm{t), following the same method, and the result is 



Sm(T + 1) 
1 ^ 

^ ' n—1 



< 



M 



(80) 



VII. Blind estimation of the spectral support 

If the spectral support of the multiband signal is unknown, 
then it is not possible to determine the set I^w, and thus the 
procedure presented in the previous sections is not applicable 
in principle. In the literature [8], [13], this problem has 
been addressed using compressed sensing techniques. For the 
SMRS scheme in Sees. |IV] and [V] this approach would lead 
to a hnear system like that in dSTl i, but in which the norm 
of the set of coefficients d,, must be minimized 



{5p} — argmin \5p\ subject to A^^, 



P=Pi 



P2 

p=Pi 



(81) 



Here, {Sp} denotes the set of coefficients 6p, and the frequency 
indices are assumed to lie in [Pi, ^2]- However, the multiband 
signal z{t) can be approximated in several intervals R{Th,T), 
h = 1,2, . . . , H, with distinct Th using polynomials Zu,(i; r/j). 
These polynomials have the same index support Izw and, 
besides, the can be selected in the infinite SMRS scheme, 
so that the relative sampling instants tk^q are the same for 
all of them. In this setting, the blind estimation problem can 
be viewed as a compressed sensing problem with multiple 
measurements vectors (MMV), [24]. Nevertheless, since the 
duration of the multiband signal can be arbitrarily large, the 
blind problem can be posed using any of the techniques from 
direction of arrival (DOA) estimation, [25, chapters 8 and 9]. 
This alternative solution is analyzed in the next sub-section, 
where the MUSIC algorithm is adopted. It is then tested 
numerically in Sec. IVIII-CI 



A final aspect is whether the SMRS scheme leads to a full 
column rank linear system, for any set I^^. In the literature, a 
sampling scheme of this kind is termed "universal", [8], [13]. 
It turns out that the finite SMRS scheme in Sec. |IV] can be 
converted into a universal one in a simple way. To see this, 
recall the explicit expression for Ziu{t; r) in (l35T l, 

zUt; r) « ~z^{t; r) = X Cp{T)e^^^P'/^ , (82) 

and assume that the sampling scheme is multi-coset sampling 
with repetition period T and minimum spacing At. This means 
that Zuj{t; t) is sampled at instants to + UrAt, r — 1,2, . . . , R, 
where the rir are distinct integers, RAt = T, and R > \Izw\- 
Evaluating (|82] | at the sampling instants yields 



z^{to + nrAt;T) = Cp(T)e 



]27vpto/T j2npn^/R 



(83) 



P&Izu 



This linear system may not have full column rank. However, 
note that ( |82] | is also valid if T is replaced with a T' following 
T' > T, given the norm of z^{t;T) in R{0,T') is 
larger than in R{0, T), and the interpolation accuracy was 
already enough in R{0, T). To replace T with T' may slightly 
increase the number of elements of \Izw\, i-e, the spectral 
sampling will be finer, and in general there will be two 
different sets Izw{T) and Izw{T'). However, if the difference 
R — \Izw{T)\ > is large enough, then one may expect that 
it is also R — \Izw{T')\ > 0, i.e, the linear system for the 
larger period T' will also have more equations than unknowns. 
Knowing this, take T' = AtP where P is the smallest prime 
number following P > R, and evaluate ( |82] | with T' in place 
of T, 



(84) 



z,„(io+n,Ai;T)« ^ Cp{T)e^^^P*'> /t' ^j2npn^ / p ^ 
P&h^ 

If R— \Izw{T')\ > this system has full column rank, since 
all minors of trigonometric Vandermonde matrices of prime 
order have full rank, by Chebotarev's theorem; (see [26, page 
25]). So, for Zu,(f;r) one can turn a multi-coset sampling 
scheme into a universal one by slightly increasing T. The 
only condition for this is that there must be a minimum over- 
sampling, in the sense that the difference R — \Izw{T')\ must 
not be too small. 

Finally, note that the SMRS scheme in Sec. |IV] is a multi- 
coset scheme if one sets to = — t and 



At = T/iulcuikQk 



(85) 



where 1cm denotes the least common multiple and v is any 
positive integer. So, the method of perturbing T also works 
for SMRS sampling. In practice, the difference between using 
T or T' can be immaterial. For example, in Sec. IVIII-BI the 
SMR scheme uses the Q values in dlOll ) with 



K 



Y[ Qk = 159049016335440. 



(86) 



k=l 



This product is one of the possible denominators in ( I85I I for a 
specific V, and the smallest prime P following P > Hs^i Qk 



IS 



P = 159049016335453. 



(87) 
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So, the relative perturbation of T is 

T' -T P 



1 = 8.17- 10"^^ 



(88) 



T R 

which is well below the numerical precision in that section 



A. Blind sampling of a collection of sparse trigonometric 
polynomials with common index support 

Consider a collection of trigonometric polynomials 
ai(t; J), a2(i; J), . . . , anit; J) with common index support 
J and coefficients /S^.p of the form in (fT9] l. 



(89) 



Assume that each polynomial is sampled at the distinct epochs 
ti,t2, ■ ■ ■ ,tN- If the samples are denoted by ah,n. 



OLh,n — dh^n'i J) + lJ-h,n, 

where fih.n is an unknown perturbation, then Eq. 
the following model for ah,n, 



ah-. 



(90) 
yields 

(91) 



where p{r) is an ordering of J, i.e, p{r) runs through all 
elements of J for r = l,2,...,|J|. Next, it is straight-forward 
to convert this equation into a matrix model, similar to those 
in DOA estimation. For this, define the matrices 

[p]r ^ p{r), [0(p)]„ ^ e^-2-f*"/^, [*(p)].,, ^ 0([p],). 

(92) 

(93) 



The model is 



A = *(p)B + N. 



The purpose now is to estimate the elements of p, which 
specify the spectral support of the polynomials z^(t; Th), and 
the matrix of coefficients B. In principle, it is possible to 
apply to (|93] l a wide variety of algorithms for DOA estimation 
in array processing, [25, Chapters 8 and 9]. However, this 
estimation/detection problem has a few differences relative to 
the typical DOA problems: the length of p can be large (on 
the hundreds), the components of this vector must be integers 
and, as will be shown in the numerical examples, the problem 
of detecting the number of components of p is not critical. 

An algorithm that seems well suited for ( |93] l is MUSIC, 
[27]. If P denotes an upper bound for \Iziu\, and U is an 
N X [N — P) matrix that spans the noise subspace of A, 
U^U = I, then the normalized MUSIC spectrum is 



x{p) 



mp)\? 



\v^<P{pW 



(94) 



The effectiveness of this estimator will be demonstrated in the 
numerical example in Sec. IVIII-CI 



VIII. Numerical examples 

A. Validation of the interpolation method 

In Sees. and |III1 a basic argument was that a band- 
limited signal which is concentrated in a time interval can 
be well approximated by a trigonometric polynomial, and this 
concentration was achieved by means of a window function. 
To test this property numerically, a BPSK signal s{t) was 
generated, whose modulating pulse was a raised cosine with 
roll-off factor 0.8, peak amplitude As — 1, and bandwidth 
B = 136/T, (chip period = 0.01324T). The formula for 
this signal is 



^ apgrc{t-pTc), 

P— — 00 



where grc{t) is the raised-cosine pulse 

. / t X cos{nl3t/T,) 

S,rc{t) = TT Sine ; -77, 



(95) 



(96) 



and the coefficients Op took two real values ±A' at random. 
A' was selected so that the peak amplitude of the signal is 
As. Then, the window in Sec. IIII-Al was applied to this signal 
with = 13.6/T and S = 0.0103 at r = 0. This gives 
e = 10~^. Fig. [1] shows the BPSK signal, the window w{t), 
and the windowed signal (bold line). The BPSK signal can be 
regarded as a multiband signal with a single component and 
spectral support [—B/2,B/2]. Therefore, the results in Sec. 
Hn] apply, and there is a set of coefficients Cp(0) such that 



s(t) 



1 



PB 
P=-PB 



Cp{0)e 



j27rpt/T 



(97) 



where ps ^ [{B + Bu,)T/2\. In order to check whether 
this is an accurate model, the signal s(t)w(i) was sampled 
with rate l/(4(i? + Byj)) (over-sampling factor 4), and then 
a polynomial like that on the right of (l97T i was fitted to these 
data using least squares. Then, the resulting polynomial, but 
divided by w(t), was used to approximate the initial BPSK 
signal. Fig. |2] presents the error of this procedure. Note that 
it increases toward the interval limits due to w(t), but the 
interpolation error is very small. Actually, for t in roughly 
R{Q,T/2), it is below -200 dB. 

B. Performance in the presence of noise 

In order to test the performance in the presence of noise, 
a multiband signal with five components was generated. Its 
components were signals of one of the following types: 

1) QPSK signal modulated by a raised cosine pulse with 
roll-off factor 0.8. This signal has the form in Eq. ( |95l l, 
but with symbols in the constellation ±A'±jA'. A' was 
selected so that the peak absolute value of the signal was 
one. 

2) Sum of 150 undamped exponentials with random phases 
and frequencies. 



150 

p=i 



pj(27r/pt+0p) 



(98) 
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Fig. 4. Supports of the multiband signal before and after windowing in Sec. 
IVIII-BI The amphtude along the y-axis has no meaning. For the numerical 
values that specify these bands, see table in Sec. IVIIl-BI 

The phases were taken at random in the interval [0, 2tt], 
and the amplitudes Up in the interval [0, 1]. These 
amplitudes were then scaled so that the peak amplitude 
of the signal was equal to one. The frequencies were 
taken at random in the interval [~B /2, B /2], where B 
is the bandwidth assigned to the signal. 
3) Sum of 200 delayed sine pulses with random delays and 
amplitudes 

200 

^apSinc(B(t-Tp)). (99) 
p=i 

The delays were taken at random in R{0, 2T), knowing 
that the sampling interval is i?(0, T). The amplitudes ap 
had random absolute value in [0, 1] and random phase in 
[0, 2tt]. They were then scaled so that the peak amplitude 
of the signal was one. 
In this example the time and frequency variables were normal- 
ized so that T — 1. For the five components, the component 
index (m), type of signal (type), central frequency (/c), band- 
width (B), bandwidth after windowing (B'), and the first (/„) 
and last (/b) frequencies in the trigonometric approximation 
were the following: 



m 


type 


fc 


B 


B' 


fa 


fb 


1 


1 


308.892 


60.4428 


69.5628 


275 


343 


2 


2 


596.276 


41.7585 


50.8785 


571 


621 


3 


3 


920.824 


39.9765 


49.0965 


897 


945 


4 


1 


1169.11 


66.6665 


75.7865 


1132 


1207 


5 


2 


1381.22 


19.1557 


28.2757 


1368 


1395 



These bands are schematically plotted in Fig. H) A set of 
moduli Q was selected using the method in Sec. |IV] [around 
Eq. (l57]il. The set was 



Qk^67 + k, k^l,..., K, within = 4. (100) 

This set ensured that the linear system in ( fSTl l had full 
column rank and was minimally over-determined. Actually, 
the number of samples was above the number of unknown 
coefficients just by one; (274 samples were taken, but there 
were 273 unknowns). However, the noise factor in (|56l l was 
l{Izw) — 48.75 dB, that is, there were interpolation instants 



Fig. 5. Noise factor in example of Sec. IVllPB] 

in which the signal-to-noise ratio would be degraded 48.75 
dB, (assuming complex white noise). In order to reduce this 
figure, several sampling grids were added. The final set of Q 
values was 

11, 18, 19, 37, 49, 68, 69, 70, 71. (101) 

These additional Q values were selected so as to reduce ^{Izw) 
as much as possible sequentially. With the set in ( IIOII ), the 
noise figure was ^{Izw) ~ 18.77 dB, i.e, in the worst case the 
noise would be amplified 18.77 dB. Fig. |5] shows 'j{t,Izw)- 
Note that except for a few peaks, 7(i, Izw) is well below this 
figure. The performance in terms of bandwidths and rates was 
the following: 

Nyquist bandwidth: 1112.13 
Landau limit: 228 
Landau limit after windowing: 273.6 

Sampling rate: 394 
(Nyquist rate)/(Sampling rate): 2.82 
(Sampling rate)/(Landau limit): 1.73 

The sampling rate equaled 394, that is, it was only necessary 
to take this number of samples in each period T. The ratio 
between the sampling rate and the Landau limit shows how 
far this implementation was from the minimal sampling rate 
(factor 1.73). This extra over-sampling was produced by the 
windowing and by the over-determination of the linear system 
so as to reduce ^{Izw)- Relative to the usual Nyquist rate, this 
implementation afforded a sampling rate reduction of factor 
2.82. 

As commented in Sec. |IV] the linear system that must be 
solved for obtaining the coefficients 5p is that in dSTl ), which 
is a sparse system. In this example, the matrix coiTesponding 
to ( l57b has only 2.18% of non-zero elements which are equal 
to one. Thus, solving this system using an efficient method for 
sparse linear systems, like that in [28], could be more efficient 
than directly multiplying by the pseudo-inverse of the matrix 
corresponding to ( |57] ). 

Fig. |6] shows the module of the coefficients 5p in this 
example, assuming that the ratio between the power of the 
multiband signal and the noise was 70 dB. The noise compo- 
nent in these coefficients is also plotted in this figure. Note that 



12 



Noisy 



Noise 




500 1000 
Normalized frequency fT 



1500 



Fig. 
Sec. 



6. Noisy 
IVIII-BI 



trigonometric coefficients and their noise component in 



Error 




. 5 



Normalized time t/T 



Fig. 7. Interpolation error for the first multib and component in example in 
Sec IVIII-"B] and conesponding eiTor bound in )80) . 



the coefficients 5p are correctly estimated with signal-to-noise 
ratio equal to 54.6 dB. 

Finally, Fig. [T] shows the error in the interpolation of the 
first multiband component, together with the error bound in 
dSOl l. Note that the error is small in a wide range and that the 
bound is very conservative. 

C. Blind sampling 

In order to test the MUSIC algorithm in Sec. IVII-AI the 
example in the previous section was repeated but for QPSK 
signals (in order to reduce the computational load). Also, the 
set of moduli Qk was changed to 



20,50,70,110,130, 



(102) 



so that all moduli Qk have common factor 10. Due to this 
property, the sampling scheme had period T/10, which made 
it possible to use the central instants = Th/ 10, i.e, the 
multiband signal was windowed in H — 500 consecutive 
intervals R{Th/10,T) with 90% overlap. Fig. [8(a)] shows the 
MUSIC spectrum in Eq. ( |94l i assuming a signal subspace of 
dimension P = \Izw\- Note that this estimator is able to detect 
the multiband components as well as their approximate width. 
(Compare this figure with Fig.|4]) Figs. |8(b)| and |8(c)] show the 



same spectrum but assuming a signal subspace of dimension 
\Izw I — 50 and \Izw \ + 50, respectively. Note that the detection 
capability remains unchanged, relative to the case P = \Izw\ 
in Fig. M 

IX. Conclusions 

A method for sampling bounded multiband signals has been 
presented, that makes it possible to approach the Landau 
limit, while keeping at the same time the noise sensitivity 
at a low level. The method is based on approximating the 
product of the signal with a window function by means of 
a trigonometric polynomial. This polynomial "inherits" the 
multiband property from the signal. The associated sampling 
scheme is the recently proposed synchronous multi-rate sam- 
pling. Also, it was shown that the blind sampling of a bounded 
multiband signal can be reduced to the sampling of a collection 
of multiband trigonometric polynomials with unknown index 
support. This fact made it possible to apply a DOA estimation 
algorithm (MUSIC) for the detection of the spectral support. 
The performance of the methods in the paper have been 
validated in several numerical examples. 
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Appendix I 

Bound on the tails of the window function 

The computation of the error bound e in Eq. ( fTSl ) depends 
on bounding the series 

1 



(103) 



to{a + n)^(a + n)2-62 

for constants a > 6 > 0. For this, introduce first the variable 
change 

\a + n/ 

to obtain 

1 



(104) 



1 x(n) 



(105) 



^ y/l-x{n) 



The function — x is increasing in [0, 1[ and the maxi- 

mum value of a: = x{n) in ( 1104b is {h/af'. So, this function 
can be bounded as 

^ < ^ = " ■ (106) 



^l-x{n) - v/1 - {b/aY Va 



52 



Coming back to (I105l l, it follows that 



S{a,b)< 



1 



CO C30 



n=0 



- 62 (a + n)2 

^ 71—0 ^ ' 



(107) 
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where ip'{a) is the derivative of the polygamma function. 

Let us proceed with the derivation of the bound. RecalHng 
the definitions in Eqs. ( l38T l and ( [39] l. since the sine function 
is bounded by one, it is 

\w(t)\ < — , 

^ ^' - nSB^t 7r(l - 5)B^y/f^ - p^T^i 



1 



sinc(i(l - S)pB^T/2) 



sinc(j(l - 5)pB^T/2)n^6{l - S)Blt^t^ - p^T^A 
C- ^ 



(108) 



t^t^-p^Tyr 

where C is the constant 

C = (siiic(j(l - 5)pB^T/2)7:^S{l - S)BI)~\ (109) 
Now, 

1 

^2 {t +pT)^{t + pTf-p^TyA 



Y,\^{t+pT)\<CY,. 

p=2 



h ( + p) V{t/T+pr-py4 

c_ 

rp2 



S{2 + t/T,p/2). 



(110) 



Also, 



5^ |w(t -PT)\^J2 +P^)\^ ^S(2 - t/T, p/2). 

p=2 p=2 

(111) 

So, combining the last two inequalities, it is 



and 



|p|>i 

< ^(S(2 - t/T, p/2) + S(2 + t/T, p/2)). 



E|w(t + pT)|<|w(t-T)| + |w(t + T)| 
+ ^(S(2 - t/T, p/2) + S(2 + t/T, p/2)). 



(112) 



(113) 



The final bound is obtained by substituting ( 1107b and ( |109t 
into this inequality. 
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Periodic repetition with period T, Eq. {5) 

I J| is tlie number of elements of the finite set J 

Data matrix in )92t 

Spectral support of Sm(t) 

norm of set of samples r](tn) in j75t 
Amplitude of s(t), |s(t)l < As 
Amplitude of Smj t) 
Data value in j90t 

Generic sparse trigonometric polynoinial in )20t 

Trigonometric polynomial in Eq. j89t 

Two-sided bandwidth of s(t) 

Two-sided bandwidth of w(i). Sec. Ull 

Bernstein space of type ttB 

Coefficient of a{t\ J) in l|20) 

Coefficient of J) in Hg) 

Coefficient matrix in )92t 

Coefficient of (t; t) 

Window parameter in Eq. j38t 

Unknowns in sparse linear system, Eq. (43) 

Lower bound on w(t) in i?(0, Ti), 5„ > 0, Eq. E) 

Pseudo-inverse of the linear system in j2U 

Vector of exponentials in (92) 

Matrix of exponentials in (92) 

Supremum of sensitivity measure in J?(0, T), Eq. (56) 
Sensitivity measure in Eq. (55) or (76) 

Generic inteipolation functions in (5) 

Index set of polynomial Sw^m {t;T), Eq. GD 

Index set of polynomial Zm{t; r), Eq. )15t 

Identity matrix in Sec. IVII-AI 

Arbitrary finite sets of integers (indices) 

Pseudo-inverse of linear system in (49) 

Number of multiband components 

Noise matrix in (92) 

Integer shift in Eq. )6U 

Ordering of set J in j9U 

Vector of indices in (92) 

Time interval defined i n (4 ) 

Window parameter in (39) 

Generic band-limited signal 

Component of z(t) 

Windowed multiband component in 13 It 
Periodized version of s^j, m{t',T). Eq. (34) 
Spectral support of z(t) in (s) 

Spectral support of Zu,{t; r) = z{r + t)w(t) in Eq. )13t 

Generic sampling instants in R{0, T) 

Length of the interpolation interval, Sec. |ll] 

Length of interval R{t,Ti) in which the interpolation 

accurate. Sec. |ll] 

Instant around which z(t) is interpolated. Sec. Ull 

Matrix spanning the noi se s ubspace of A, Sec. IVII-AI 

Noise sample in model (90) 

Interpolation polynomials in j23t 

Interpolation polynomials in (53) 

Band-limited window. Sec. |ll] 

Multiband signal with bounded components 

Windowed multiband signal in (29) 

Periodized version of z»(t;T), Eq. (35) 
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